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We analyze the evolution of the Gaussian discord between two resonant harmonic oscillators 
coupled to a common environment. For this, we use the same tools we applied before to fully 
characterize the evolution of the entanglement in this system (J. P. Paz and A. Roncaglia, Phys. 
Rev. Lett. 100(2008)). The asymptotic value of Gaussian discord is obtained as a function of 
parameters characterizing the environment (temperature, couplings, etc) and the initial state of 
the system (initial squeezing, initial purity, etc). The type of Gaussian measurement optimizing 
the extraction of information between the oscillators is fully characterized by means of a phase 
diagram. Such diagram (with phases corresponding to homodyne or heterodyne measurements) has 
similar topology to the one describing dynamical phases for the entanglement. We present evidence 
pointing to the fact that Gaussian discord is not always a good approximation of true discord as 
the asymptotic value of the former is shown to be a non-decreasing function of temperature (in the 
high temperature regime), reaching an asymptotic value of log(2) for a pure initial state (and lower 
values for mixed initial states). 
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I. INTRODUCTION 

Quantum discord is a measure of quantum correla- 
tions that includes not only those contained in entangled 
states, but also others that are present in separable states 
[l[ . The fact that separable states can have non-classical 
correlations, which are quantified by quantum discord, 
came as a rather bit surprise and has raised considerable 
attention in recent years. Such attention was partly fu- 
eled by recent observations suggesting that discord may 
also be the computational resource behind the power of 
the DQC1 model of computation |2|, [3J. Quantum dis- 
cord has also been related to a variety of problems, such 
as the complete positivity of quantum dynamics Q and 
the state merging protocol^ [||, among others. The evo- 
lution of quantum discord for quantum open systems has 
also been studied for systems of qubits[7H9(, in both the 
Markovian and non-Markovian regimes. For continuous 
variable systems much less is known. This is because 
computing the quantum discord is rather hard since it 
involves solving a optimization problem that is highly 
nontrivial even for the simplest case of two qubits. How- 
ever, for a system of two particles (two modes) that are 
restricted to Gaussian states it is possible to compute an 
approximation to quantum discord. This is the so-called 
Gaussian discord, which has been defined and computed 

recently HH EH . 

Here, we fully characterize the behavior of Gaussian 
discord for a system of two harmonic oscillators coupled 
to a common environment. The environment is itself a 
collection of independent harmonic oscillators (linearly 
coupled with the system). This same model was used 
before to study entanglement dynamics [Pa. Il3| . where 
a complete characterization of the entanglement in the 
equilibrium state was given. Our analysis is exact and 
includes both non-Markovian and non-perturbative ef- 



fects. In particular, we analyze here the same scenario 
studied before by Maniscalco et al. in Ref. [14j . where 
the evolution of Gaussian discord was analyzed using a 
weak coupling approximation (that, as we see, prevents 
one from observing rather interesting behavior such as 
the existence of entanglement at long times for initially 
uncorrelated states). Also, we do not restrict our analy- 
sis (as was done in Ref. [Tj]) to high temperatures and 
include non-Markovian effects well beyond the short time 
regime. 

The paper is organized as follows: In Sec. |H] we de- 
scribe the model and the technique, based on the use 
of the exact master equations. In Sec. IIIII we discuss, 
for completeness, quantum discord and its Gaussian ap- 
proximation. In Sec. II VI we present our main results. We 
include a full characterization of the behavior of Gaus- 
sian discord in the asymptotic state and study features of 
the temporal evolution. We also compare the dynamics 
of Gaussian discord and entanglement. Finally, Sec. IVl 
contains the conclusions of the work. 



II. THE MODEL 

We consider a system of two harmonic oscillators 
with coordinates (xi,pi) and (£2,^2) with the following 
Hamiltonian (that includes both position and momentum 
coupling): 
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This system interacts with an environment formed by 
N independent harmonic oscillators (with a Hamiltonian 
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environment interaction is bilineal: 
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The total Hamiltonian, H = H sys + H env + Hi ntl is 
quadratic in position and momentum. Then, tempo- 
ral evolution maps initial Gaussian states onto Gaussian 
states. Throughout the paper we will consider initial 
states po with no correlations between the system and 
the environment: pa — ps ® Pe- Moreover, we will as- 
sume that the initial state of the environment pe is of 
thermal equilibrium at certain temperature T. The effect 
of the environment on the system is fully characterized 
by the initial temperature T and by the spectral density, 

defined as J(ui) — J2 n =i ^( w — u n) 2 Cn u • We cons ider a 
family of spectral densities: 
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where 0(x) is the step function. In the continuum limit 
the environment contains oscillators with frequencies up 
to the cut-off A. The environment is Ohmic for n = 1, 
sub-Ohmic for n < 1 and super-Ohmic for n > 1. The 
strength of the coupling is controlled by the parameter 

7o- 

As in Ref. [l3|, we consider a family of system- 
environment interactions with a certain symmetry. This 
is evident if we rewrite the Hamiltonian in terms of nor- 
mal modes, whose coordinates are related to those of 
the original modes as: x± — -75(^1 ± #2) and p± — 
-75 (pi ± Vi)- Due to the symmetry of the coupling, the 
interaction with the environment involves only the x+ 
and p+ operators. In terms of these modes the Hamilto- 
nian of the system is: 
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where m± = 

cj 2 (1±^j|-) (l ± fjf- J are the masses and frequen- 
cies of the new modes, respectively (notice that + and 
— modes are normal modes since the two oscillators 
are assumed to be resonant). Two relevant cases can 
be distinguished within the family of couplings we 
considered: (1) When the coupling is only trough the 
position coordinates, i.e., c\i = c n = 0. (2) When 
the coupling is symmetric in position and momentum: 
c n = c n . The master equation for the state of the system 
arising for each coupling is well known and is reviewed 
for completeness in the next subsection. 



A. Position coupling 

In this case, the evolution of the reduced density ma- 
trix of the system, p, is ruled by the exact master equa- 



tion for Quantum Brownian Motion [T^, LL5| : 

p= -i[H R ,p] -ij(t)[x + ,{p + ,p}]- 
-D{t)[x+,[x+,p]]-f{t)[x,]p+,p]] 



(5) 



Where Hr is the renormalized Hamiltonian of the sys- 
tem (H R = H s + fSfl 2 (t)x 2 + ). The action of the en- 
vironment is fully contained in the time dependent co- 
efficients 5ft 2 (t), j(t), D(t) and /(£). In turn, these 
coefficients can be expressed as functions of the initial 
temperature T of the environment and its spectral den- 
sity. The environment induces a renormalization on the 
frequency of the oscillators. The frequency shift 5fl (t) 
has a long time limit that depends on the spectral den- 
sity only (and depends on the cut-off frequency A). In 
what follows we choose the physical parameters in such a 
way that the long time properties are cut-off indepen- 
dent. For this, we must choose 'bare' parameters u> 2 
and C12 accordingly (to absorb the cut-off dependency 
carried by the frequency shift induced by the coupling). 
For example, in the case of two non-interacting oscil- 
lators one must take C12 = — lim t _ i . 00 5Vt 2 (t) so that 
limt->oo C12 = 0. For the spectral density ([3|) one can 
show that limt->oo ^ 2 (*) = — - ^77 ■ The equilibrium 
value of the coefficients will be noted simply as f2, C\%, 
7, D, etc. The temporal dependence will be written ex- 
plicitly if needed. 

The Gaussian nature of quantum states is preserved 
by time evolution. This simplifies the description of 
the evolution since the full state is specified by the 
first and second moments of the quadrature operators 
R = {xi,pi,X2,p2), i-e., by their mean values (Hi) and 
the covariance matrix &ij = ^({Ri,Rj}) — (Ri)(Rj). 
From the exact master equation ([5]) it is easy to obtain 
the following evolution equations: 
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Here, Q+(t) is the renormalized frequency of the + mode. 
The evolution equations for the — mode are simply those 
of a free harmonic oscillator, and describe a rotation 
in the phase space by an angle ui—t. From the previ- 
ous equations it is clear which is the role of the time 
dependent coefficients of the master equation: j(t) is 
a damping coefficient that induces loss of energy, D(t) 
is a diffusion coefficient that increases the dispersion in 
momentum (and therefore also in position), and fit) is 
an anomalous diffusion coefficient that affects only the 





dispersion in position. This last coefficient reflects the 
asymmetry of the interaction between system and envi- 
ronment. 

The asymptotic state is attained after a timescale fixed 
by the long time value of j(t) (of course, this happens 
after the time dependent coefficients approach constant 
values, which is the case for typical spectral densities such 
as the ones we consider here). Thus, from (|7|) one can 
obtain the asymptotic values of the dispersions Ax + = 

2 + ) and Ap 



mf2 



Also, the correlations between x+ and p+ vanish in 
the asymptotic limit: ■^{{ X +;V+}) = 0. From Eqs. [5] 
it is clear that due to the anomalous diffusion term 
the equipartition principle is violated (or, equivalently, 
Ap + ^ mfl + Ax+), i.e., there is squeezing in the asymp- 
totic state. This fact have direct consequences in the 
asymptotic behavior of the correlations between the os- 
cillators. 



B. Symmetric coupling 

When the system-environment interaction is symmet- 
ric in position and momentum (c„ = c„), the master 
equation is [12| 

f>(t) =-t[Hr,p] -il(t)([x+,{p+,p}] - \p + ,{x+,p}})- 



-D(t)([x + ,[x+,p}} + 



t\p+,\p+,p1 
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Now, the renormalized Hamiltonian is: 



Hfi — Hs 



sn 2 (t) 



p-\ 



m+x 2 . 



2m+(jji 



Again, the time dependent coefficients depend on the 
spectral density (D(i) also depends on the tempera- 
ture). As expected, the operators x+ and p+ appear 
in a symmetric way in the master equation (and, as a 
consequence, there is no anomalous diffusion term). The 
asymptotic dispersions are: 



m£l+Ax+ — Ap + = \ — 



(10) 



Where D and 7 are the asymptotic values of D(t) and 
7(t). Also, ({x+,p+}) = 0. There is no squeezing in the 
asymptotic state of the + mode (i.e., the equipartition 
principle is always fulfilled). 



The renormalization in the parameters of the system 
is as follows: 
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We considered that the modes 1 and 2 do not interact in 
the long time limit (C12 = C12 — 0)- With that choice 
M = to_ and fi = w_ . 



C. Asymptotic state 

The master equation enables us to obtain not only the 
asymptotic covariance matrix of the + mode. In fact, in 
the same way we can show that all correlations between 
the ± modes vanish and that the — mode evolves freely. 
It's covariance matrix, in general, will depend on time. 
Therefore, the covariance matrix corresponding to the 
asymptotic state has the following form in terms of the 
± coordinates: 
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The first 2x2 block on the diagonal is the covariance 
matrix of the + mode. The second 2x2 block on the 
diagonal is the covariance matrix of the — mode, which 
is obtained as the transformation of an initial diagonal 



covariance matrix cr 
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under the 



rotation in phase space corresponding to the free evolu- 
tion: 



-(t) = W)<r-(0)E T (t) 
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From this, we can obtain the covariance matrix for the 

original modes. It has the form er 12 = ( t ~ J with 

a = j3, reflecting the symmetry of the model. The com- 
ponents of ar 12 have simple expressions in terms of Ax+, 
Ap + , Ax-, Ap_ and t, but they are not all necessary. In 
fact, any measure of correlations must be invariant un- 
der local unitary operations. If these operations are also 
Gaussian, they correspond to local symplectic operations 



at the phase space level. Therefore, the only quantities 
one need to know to evaluate measures of correlations 
like the mutual information or logarithmic negativity are 
the local symplectic invariants of the covariance matrix. 
These invariants are simply the determinants of each 2x2 
block and the determinant of the entire covariance ma- 
trix. The local symplectic invariants of er 12 admit this 
compact form: 



their nature. Classical correlations are quantified by the 
maximum amount of information J(A : B) that can be 
obtained about a subsystem, say A, by means of local 
measurements performed on the other: 
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B = dct(/3) = dct(a) = A 
C = dct( 7 ) = 1 (4> 2 + 
D = det(tr) = 4>\4> 2 _ 



where 



Ax+Ap + and 



purity of each mode. All the temporal dependence is 
contained in the function h(t,r,r cr i t ): 



where {II m } is a positive operator-valued measure 
(POVM) corresponding to a measurement in subsys- 
tem B. The probability of obtaining the result m is 
(15) P" 1 = ^ r (/°sn m ). The state of A after obtaining this 
result is PA\m = TrB(pAB^m)- Then, the sum in the 
last expression is the average entropy of subsystem A af- 
ter a measurement described by the POVM {II m } is per- 
formed on B. Classically, this sum would coincide with 
Ax_Ap_ measure the the conditional entropy S(A\B) = S(A,B) — S(B) and 



+ <f> 2 _ - 2cj) + (j)-h(t, r, r crit )) 



*-, 



h(t, r, rcrit) = cos 2 (w_i) cosh(2(r - r mt ))+ 
sin 2 (w_£) cosh(2(r + r cr it)) 

where r and r cr i t are the squeezing of each mode: 
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From Eq. (TO)]) it is evident that if the squeezing of either 
mode is zero, then h is constant and there is no temporal 
dependence in the symplectic invariants. 

Summarizing, any measure of correlations in the 
asymptotic state can be evaluated from the four local 
invariants A, B, C and D of Eq. (IT5|) . The only mem- 
ory about the initial state is contained in the — mode, 
through </>_ and r. On the other hand, the + mode 
reaches the thermal equilibrium with the environment, 
so <p + and r cr n depends exclusively on the temperature 
T and the characteristics of the environment and its cou- 
pling with the system. 



III. QUANTUM DISCORD 

Quantum discord arises as the difference between two 
classically equivalent expressions for the mutual infor- 
mation of a bipartite system. It can also be interpreted 
as the difference between total and classical correlations. 
Therefore it is a measure of quantum correlations. To- 
tal correlations contained in a given state pab of a bi- 
partite system are quantified by the mutual information 
I(A :B): 



I(A : B) = S(p A ) + S(p B ) - S(pab) 



(18) 



where S(-) is the Von-Neumman entropy and pa(b) is the 
reduced state of the subsystem A(B). Mutual informa- 
tion takes into account all the correlations, regardless of 



then I(A : B) would be equal to J (A : B) . The quan- 
tum discord D(A : B) is the difference between these 
two quantities: 



*5(A : B) = I(A 



B) - *J{A : B) 



= S{p B ) - S(pab) + mf I ^p m S{p A \ m ) J 

(20) 

The main difficulty in calculating the quantum discord 
of a given state is the optimization over all the possibles 
measurements in one subsystem. So far, for the simplest 
case of two qubits, an analytical expression for quantum 
discord is only available for a reduced family of states. 
For continuous variable systems there is only an approxi- 
mation to the quantum discord, called Gaussian quantum 
discord, developed by Adesso and Datta[l0| and simul- 
taneously by Giorda and Paris |lll| in 2010. The main 
restriction of this approximation is that the optimiza- 
tion over all the measurements over one single mode is 
restricted to the set of Gaussian measurements, or gen- 
eralized Gaussian POVM's. For this reason the Gaus- 
sian quantum discord provides only an upper bound for 
the unrestricted quantum discord (except for a family of 
states, for which Gaussian measurements are known to 
be optimal). In the following we briefly review the main 
points needed to evaluate the Gaussian quantum discord 
of a two mode Gaussian state. 



Two mode Gaussian states 



In what follows we consider dimensionless quadrature 
operators R — (x\,pi,X2,P2), i-e., Xi — -js{a.i + a\) and 
-^| (a,i + a\ ) where a^ is the annihilation operator of 
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each mode. The Gaussian approximation to the quantum 
discord applies to two mode Gaussian states. These are 
completely described (within local displacements) by the 
4x4 covariance matrix er 12 of the quadrature operators 



°"12 — 



a 7 
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where a and (3 are the 2x2 covariance matrices for the 
modes 1 and 2, respectively, and 7 is the matrix contain- 
ing the correlations between (xi,p{) and (X2,P2)- As was 
noted in Sec. Ill CI it is necessary only to know the local 
symplectic invariants of the covariance matrix. They are 
A = deta, B = det/3, C = det7 and D = dot er 12 . Fur- 
thermore, by means of a symplectic operation (generally 
global), any Gaussian state can be transformed in an- 
other Gaussian state with covariance matrix of the form 
a = diag{v +1 v + , i/_, f_). The symplectic eigenvalues v+ 
and V- are invariants under symplectic transformations 
and can be calculated as 2v 2 _ — A±\/A 2 — AD, with A = 
A + B + 2C. A given covariance matrix represents a valid 
physical state if and only if v± > 1/2, that is, if the uncer- 
tainty principle is satisfied. Finally, the entropy of a state 
with symplectic eigenvalues {v{\ is S — Y^iffai)) with 
/(a?) = (x + 1/2) log (as + 1/2) + (x - 1/2) log (x - 1/2). 



B. Gaussian quantum discord 



The only non trivial point in evaluating the expres- 
sion (12TJ1) for a two-mode Gaussian state is the optimiza- 
tion of the last term over all the possible measurement 
on the mode 2. The approach of Refs. [T(J, QM con ~ 
sists in restricting the optimization to the set of Gaus- 
sian POVM's. The elements of a Gaussian POVM can 
be expressed as 11(77) = ^^MPoW^y where p is a 
valid density matrix of a single mode Gaussian state, 



W{rf) 



r ' a i is the displacement operator and r\ 



is a complex number corresponding to each of the pos- 
sible results. The POVM is completely determined by 
po- If the global state is pn, the probability density 
of obtaining the result ry is p{r\) — Tr{p2^A{ii)) and 
the state of the mode 1 after obtaining that result is 
Pi\ v = Tr 2 (/9i2(l <8> II(?7))). Therefore, the generalization 
of the sum in (j2"0j) is f p(r])S(pi\ ri )d 2 r]. It can be shown 
that the covariance matrix e of pi\ v is independent of r\ 
and equals a — 7(/3 + <t )~ 1 7 T , where (Tq is the covari- 
ance matrix of the state po that defines the POVM. Then 
S(Pi\ri) — /(det(e)) and the last integral is trivial. Since 
/(•) is a growing function, the final step is to minimize 
E = dct(e) over all the covariance matrices er . The 
result is as follows: 



S(pia) = f(y/B) - f(v+) - /(»/_) + f(y/E~j (21) 
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FIG. 1: . Phase diagram for the asymptotic entanglement. 
Ohmic environment and position coupling. A = 20, 7 = 0.1, 
m = 1, fi = 1, C 12 = 0. 



Where E m i n is given by: 
Bmin = mincfet(e) = 



2C 2 + (l/4-S)(A-4D)+2|C|^C 2 + (l/4-S)(A-4D) .„ 
4(1/4-5)^ " 9 S U, 



AB-C 2 +D-y/C 4 + (-AB+D) 2 -2C 2 (AB+D) 
2B 



if g> 
(22) 



The quantity g that discriminates the two cases in the 
last Eq. is: 



g=(D- ABf - (1/4 + B)C 2 {A + AD) 



(23) 



The case g > corresponds to states for which the Gaus- 
sian discord is optimized by homodyne measurements. 
These measurements are projections onto pure states of 
infinite squeezing. For example, the measurement of po- 
sition or momentum is homodyne. In the other hand, 
g < corresponds to heterodyne measurements. These 
are projections onto squeezed thermal states. 



IV. RESULTS 

This section is divided into three parts. First we review 
previous results on entanglement dynamics. In the sec- 
ond part we analyze the Gaussian discord of the asymp- 
totic state. Finally we discuss the time evolution of dis- 
cord for intermediate times. 



A. Asymptotic entanglement 



Entanglement dynamics 
Il3j . A good entanglement 
rithmic negativity, defined 
where V- is the smallest 
covariance matrix that is 
position. This eigenvalue 



was investigated in Refs. [12], 
measure is given by the loga- 
as Em — max{0, — log(2i>_)}, 
symplectic eigenvalue of the 
obtained after partial trans- 
can be calculated as 2i>_ = 



A - ^/A 2 -4det(er 12 ), A = A + B - 2C. As shown in 
Refs. [12|, [l3|, the logarithmic negativity of the asymp- 
totic state is: 



where E{t) is: 



E N (t) = max{0,.E(£)} 



E(t) =E N + AE N G(t) 



(24) 



(25) 



The function G(t) is an oscillatory function with period 
7r/u;_. The mean value En and the amplitude AEn of 
E(t) are: 



E N = max{|r|, \r crit \} - S crit 
AE N = min{ | r |, | r crit \ } 



(26) 



The squeezing factors r and r cr n are defined in Eq. (1171) , 
and S cr it is given by S crit = \ \og{AAx + Ap + Ax_Ap_). 
From this, it is clear that we can identify three differ- 
ent 'phases' for the asymptotic entanglement. First, if 
E(t) < the asymptotic state is not entangled. This 
phase is called SD (Sudden Death). On the other hand, 
if E(t) > the asymptotic state is entangled. This phase 
is called NSD (No Sudden Death). Finally, there is an 
intermediate situation in which the function E(t) alter- 
nates between positives and negatives values. In this case 
the asymptotic entanglement suffers successive events of 
'sudden death' and 'sudden revivals', and so this phase 
is called SDR (Sudden Death and Revivals) . These three 
phases are illustrated in Fig. [T] for the case of position 
coupling with an Ohmic environment. For any initial 
squeezing r, there exists a temperature above which the 
asymptotic state enters the SD phase, i.e., it is separable. 
This critical temperature grows exponentially with the 
initial squeezing r, which is a resource to create entangle- 
ment. For low temperatures a surprising result emerges 
(as it is non-perturbative, it is missed by perturbative 
treatments |14|): There is an NSD island which is due to 
the fact that the asymptotic state is squeezed. For sym- 
metric coupling the phase diagram is much simpler: Only 
two phases (SD and NSD) exist. In what follows we will 
see that the dynamics of discord shows many differences 
and some similarities with the above results. 



B. Asymptotic Gaussian discord 

1. Symmetric coupling 

The Gaussian discord of the asymptotic state can be 
evaluated using Eqs. ([2"I]). (|2"2j). and $Z3§ together with 
the expressions for the symplectic invariants in Eq. (|15p . 
For symmetric coupling the + mode is not squeezed 
{fcrit = 0) and all the symplectic invariants become time 
independent (i.e., h(t,r,r crit ) = cosh(2r)). For simplic- 
ity, we consider the non interacting case (C12 = C12 = 0) 
and restrict ourselves to initially pure states for the mo- 
ment (0_ = 1/2). There are two interesting aspects of 
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FIG. 2: Type of measurement that minimizes the Gaussian 
discord of the asymptotic state. HM: homodyne measure- 
ments. HT: heterodyne measurements. Ohmic environment 
and symmetric coupling. A = 20, 7 = 0.1, M — 1, fi = 1, 
C12 = C12 = 0. 



quantum discord in the asymptotic regime. First, we 
analyze what measurement optimizes the extraction of 
information between the two components of the system 
(this is needed to obtain the actual value of the discord). 
After that, we look at the asymptotic value of Gaussian 
discord as a function of squeezing and temperature. 

Figure [2] shows the type of measurement, homodyne 
or heterodyne, that minimizes the Gaussian discord of 
the asymptotic state. It is clear that this phase di- 
agram has the same topology than the one character- 
izing the dynamics of the entanglement in the asymp- 
totic regime. For any given temperature, there is a de- 
gree of squeezing after which the optimal measurement 
is always homodyne (HM). Contrarily, heterodyne (HT) 
measurements are optimal for temperatures above a cer- 
tain curve T'(r). This is a natural result, in fact in the 
limit of infinite squeezing there is an observable associ- 
ated with the relative motion between the two oscillators 
that is perfectly localized. For example, if the observable 
g_ = -75(^1 — £2) is perfectly localized and the position of 
mode 2 is measured, then the state of mode 1 collapses to 
an eigenstate of x\ . This is a pure state with zero entropy. 
Therefore, the (homodyne) measurement of the position 
in mode 2 minimizes the discord since the minimum pos- 
sible value of the sum of Eq. (|2T)|) is achieved. It is worth 
pointing out that because the asymptotic state for our 
system is symmetric, the discord doesn't depend on the 
choice of the subsystem that is measured, as is generally 
the case. A simple expression for the curve T'[r) can 
be obtained from Eqs. {23]) and ([IS). With 0_ = 1/2 
and r cr it — 0, the discriminant g of eq. (|23[) is a poly- 
nomial of degree 8 in the variables </>+ = h coth(^) and 
h = cosh(2r). The curve T'(r) correspond to points for 
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FIG. 5: Type of measurement that minimizes the Gaussian 
discord for the asymptotic state. Position coupling and Ohmic 
environment. A = 20, 7 = 0.1, m = 1, Q — 1, C12 = 0. 



FIG. 3: Gaussian discord for the equilibrium state. Ohmic 
environment and symmetric coupling. A = 20, 7 = 0.1, M = 

1,0=1, C*12 = C*12 = 0. 
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FIG. 4: Gaussian discord for the equilibrium state. Ohmic 
environment and symmetric coupling. A = 20, 7 = 0.1, M = 
1,0 = 1, da = C12 = 0. 



which g = 0. This happens if and only if h — h C Ht{4'+)i 
with: 



-1 



I'critW+i 



l&(j)\ + 56 + </>+ 2 



1 



1204 



(27) 



The function T'(r) is defined implicitly by cosh(2r) = 
h C rit(4>+(T')). For high temperature or squeezing cf> + ~ 



g, /i crit ~ ^ and cosh(2r) ~ ^-. Then T'(r) ~ ^e 2r . 
It is simple to show that discord can never disappear 
in a finite time (i.e., there is no sudden death of discord). 
Figure |3] shows the asymptotic value of the Gaussian dis- 
cord as a function of the initial squeezing r and the tem- 
perature T . As expected, for low temperatures Gaussian 
discord grows with the squeezing. On the other hand 
the temperature dependence is rather surprising: There 
is a temperature after that Gaussian discord grows with 
temperature approaching a saturation value which can 
be shown to be equal to log(2). This is not expected 
as discord is a measure of quantum correlations, which 



decrease with temperature. This result can be analyti- 
cally found for high temperatures: Thus, in such a case 

the symplectic invariants are A = B ~ C ~ -j-, and 
D = <j)%<j) 2 _. Then, E min ~ 4(c/>_ + 1/4) 2 and the Gaus- 
sian discord is: 

D SUP = / fo+/2) - / (0+) - / (<£_) + / (2 (<£_ + 1/4)) 

(28) 
As / (<fi+/2) — f (4>+) —t — log(2), for initially pure states 
(for which 0_ = 1/2), we have D sup = log(2), as Fig. 0] 
indicates. This high temperature result for Gaussian dis- 
cord is independent of the asymptotic squeezing of the 
+ mode. Therefore, it is also valid for the case of po- 
sition coupling that will be discussed below. For high 
temperatures the optimal Gaussian measurement is al- 
ways heterodyne and in this case the non-orthogonality 
of Gaussian states seem to pose a constraint on the max- 
imal information that one can gather about one system 
by performing Gaussian detection on the other. It seems 
that this is the cause of the failure of Gaussian discord as 
a good approximation to true discord (which, as a mea- 
sure of quantum correlations, must decrease with growing 
temperature in the high temperature regime). 



2. Position coupling 

If the system-environment interaction is through po- 
sition the asymptotic state is squeezed (the squeezing is 
measured by r cr u). As a consequence, the local symplec- 
tic invariants of the asymptotic state become time de- 
pendent (they oscillate periodically with frequency w_, 
as seen in Eqs. (fT5"j) and (|TB)) ). The type of measurement 
that optimizes the extraction of information between the 
two subsystems changes because of this fact. This is seen 
in Figure [5] As opposed to the case discussed above, 
the phase diagram characterizing the optimal measure- 
ment shows three distinct phases. Again, large squeezing 
and high temperatures are, respectively, associated with 
HM and HT measurements. However, the fact that the 
asymptotic state is squeezed due to the coupling induces 
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FIG. 7: Type of measurement that minimizes the Gaussian 
discord for a mixed initial state with 4>- — 3/2. Symmetric 
coupling and Ohmic environment. A = 20, 7 = 0.1, M — 1, 

n = 1, C12 = C12 = 0. 



FIG. 6: Peak-to-peak amplitude for the oscillations of the 
Gaussian discord in the equilibrium state. Position coupling 
and Ohmic environment. A = 20, 7 = 0.1, m — 1, O = 1, 
C12 = 0. 



a low temperature phase where the optimal measurement 
is homodyne. The existence of an explicit time depen- 
dence in the asymptotic state reflects on the fact that 
the HM and HT phases are separated by a new phase 
where the optimal measurement depends on time (for 
some times the optimal measurement is homodyne while 
for other times such measurement is heterodyne). This 
is rather similar to what happens with the phase dia- 
gram characterizing the evolution of the entanglement in 
the long time limit. The curve T'(r) that separates the 
two phases for the case of symmetric coupling is plotted 
for reference. The new time dependent phase disappears 
for high temperatures (this is due to the fact that for 
the Ohmic environment we are considering r cr i t —> for 
increasing temperature). 

The asymptotic value of discord in this case is also 
time dependent: Gaussian discord oscillates with u;_ fre- 
quency The mean value of the oscillations behaves es- 
sentially like the Gaussian discord for the symmetric cou- 
pling case (Fig. [3]). Figure [5] shows the peak-to-peak am- 
plitude AD of the oscillations of the Gaussian discord in 
the equilibrium state. Again, because r cr u goes to zero 
rapidly for increasing temperature, the amplitude AD of 
the oscillations becomes negligible for temperatures that 
are a few times larger than $7. 



3. Mixed initial states 

All the previous results were obtained for pure initial 
states, i.e., 0_ = 1/2. Now we explore how the results 
change if this condition is relaxed. For simplicity, only 
symmetric coupling and Ohmic environments are consid- 
ered. 

The nature of the measurement minimizing the Gaus- 
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FIG. 8: Gaussian discord of the asymptotic state versus the 
temperature of the environment. Mixed initial state with 
4>- = 3/2. Symmetric coupling and Ohmic environment. 
A = 20, 7 = 0.1, M = 1, fl = 1, C12 = C12 = 0. 



sian discord is shown in Fig. [7] for a mixed initial state 
with 0_ = 3/2. There are obvious differences between 
this and the corresponding figure for a pure initial state 
(see Figure ^ . Some of the results are also counter in- 
tuitive: for example, in the limit of large squeezing the 
optimal measurement turns out to be heterodyne. The 
asymptotic value of discord decreases when increasing 
the initial mixedness. This is shown in Figure [5] where 
the asymptotic value of the Gaussian discord is shown as 
a function of the temperature for different values of the 
initial squeezing (for 0_ = 3/2). The asymptotic value 
of Gaussian discord in the high temperature limit can 
be analytically estimated using equation 1281 Thus, the 
asymptotic value of Gaussian discord is examined as a 
function of the initial mixedness of the state (which is 
estimated by the effective temperature derived from <fi- , 
To, which is such that 4>- = \ cott^^p)). This is shown 
in Figure [SI 
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FIG. 9: High temperature value of the Gaussian discord of the 
asymptotic state versus the initial temperature of the system. 
The analytical results (Eq. I|28|ll are compared to numerical 
simulations. 



FIG. 11: (Color online) Gaussian discord and entanglement 
evolution for an entangled initial state with r — 1. Symmetric 
coupling and Ohmic environment at T = 4il and T = 30S1. 
A = 20, 7 = 0.1, M = 1, fi = 1, C12 = C12 = 0. 
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FIG. 10: (Color online) Gaussian discord and entanglement 
evolution for an entangled and a separable initial state, both 
with initial squeezing r — 1. Symmetric coupling and Ohmic 
environment at T = 0. A = 20, 7 = 0.1, M = 1, fi = 1, 
C12 = C12 = 0. 



C. Dynamics of Gaussian Discord 

Temporal evolution of the Gaussian discord can be an- 
alyzed by a simple exact numerical solution of the dy- 
namical equations of the covariance matrix. Results are 
presented below for two sets of initial conditions corre- 
sponding to separable and entangled states for the origi- 
nal modes. 



1. Symmetric coupling. 

From our previous discussion it is clear that discord in 
the asymptotic state depends only on the initial squeez- 
ing \r\ of the — mode. This is illustrated in Fig. [TU1 where 
it is shown that Gaussian discord behaves in the same 
way for both initially separable and entangled states with 
the same degree of squeezing (r = 1 in this case). Three 
temporal scales are relevant for the evolution of quantum 
correlations. For example, for an initially entangled state 
decoherence leads to a rapid decay of both entanglement 



FIG. 12: Gaussian discord for two initial states with positive 
and negative squeezing. Position coupling and Ohmic envi- 
ronment at T = 0. A = 20, 7 = 0.1, to = 1, fi = 1, C12 = 0. 



and Gaussian discord. This takes place in a timescale 
that is much shorter than the dynamical time of the sys- 
tem t = 2ir/Cl. After this, quantum correlations oscillate 
with a decreasing amplitude. The timescale character- 
izing the oscillations is fixed by $7 while the amplitude 
decay is governed by 7. For the case of a T — envi- 
ronment, the evolution of Gaussian discord and entangle- 
ment is essentially the same. This is not true for higher 
temperatures, as Fig. [Til shows. As seen in the Figure, 
even though there is sudden death of entanglement, the 
Gaussian discord reaches an non-zero asymptotic value. 
For sub-Ohmic environments the results are similar, and 
the only difference is that the equilibrium state is reached 
faster. This is because for a sub-Ohmic environment the 
coupling with the low frequency oscillators in the envi- 
ronment is stronger. 



2. Position coupling. 

As discussed before, when system and environment 
interact through position, there are oscillations in the 
asymptotic state. The evolution of the Gaussian discord 
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FIG. 13: Gaussian discord evolution for an initial separable 
state with r = 1 and different temperatures for the environ- 
ment. Position coupling and Ohmic environment. A = 20, 
7 = 0.1, m= 1, fi = 1, C12 = 0. 
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FIG. 14: Gaussian discord and entanglement evolution for an 
initial state and temperature corresponding to the SDR phase 
(T = 0.35S1 and r = 0.15). Position coupling and Ohmic 
environment. A = 20, 7 = 0.1, m = 1, fl = 1, C12 = 0. 



for entangled initial states is shown in Fig. [T2j As be- 
fore, it is seen that the mean value and amplitude of the 
oscillations depends only on the initial squeezing. For 
this case, the asymmetry between position and momen- 
tum can be seen by looking at the evolution of initial 
states that are squeezed either along position or along 
momentum. In fact, as shown in the Figure, initial states 
with momentum squeezing (r < 0) are more sensitive 
to the environmental interaction during the initial time. 
Such states decohere more rapidly than states that are 
squeezed along position (r > 0). For those states to deco- 
here, dynamics has to take place (during which position 
and momentum exchange roles). This is clearly seen in 
Fig. 1121 The position-momentum asymmetry is also seen 
in the 7r/2 phase difference in the long time oscillations 
of both entanglement and discord. 

Figure [13] shows the evolution of Gaussian discord for 
initial separable states and different temperatures. As 
the temperature increases the amplitude decreases, in ac- 
cordance with the results of Sec. IIVB 21 Finally, Fig. Q3] 
compares the evolution of entanglement and Gaussian 
discord for an initial state and temperature correspond- 



ing to the SDR phase (see Fig. [T]). It is interesting to 
note that the Gaussian discord attains local maximum 
values when the entanglement disappears. 



V. CONCLUSIONS 

We analyzed the dynamics of Gaussian discord of a 
bipartite system of two harmonic oscillators that inter- 
act with a common bosonic environment. For a coupling 
that is bilinear (in position and momentum of both the 
system and the environment) the Gaussian measurement 
optimizing the extraction of information between the two 
subsystems can be obtained and used to compute the 
asymptotic discord. In this way, the evolution of Gaus- 
sian discord can be fully characterized in the long time 
limit. Finally, by means of a numerical exact solution of 
dynamical equations we analyzed the time dependence 
of discord for intermediate times (before the asymptotic 
state is reached). In this way, we confirmed analytical 
evaluation mentioned above. The results have some sur- 
prising features. Previous works on entanglement dy- 
namics showed that such quantum correlations can be 
created by the interaction with the common reservoir 
(in fact, non-entangled initial states exhibit long time 
resilient entanglement). However, this is the case only 
for temperatures that are low enough (for any value of 
the initial squeezing there is a critical temperature above 
which the final state is separable) . This is not what hap- 
pens with the type of quantum correlations characterized 
by discord. Gaussian discord can also be created by the 
interaction with a common reservoir. However, once it is 
created by the interaction with the environment, it never 
disappears. As shown in our work, the asymptotic value 
of Gaussian discord may be small but never vanishes (as 
expected). The fact that the asymptotic value for Gaus- 
sian discord does not decrease (but saturates) when the 
temperature grows is far from intuitive. Most likely, this 
is an indication of the failure of the Gaussian approxi- 
mation for discord in this case. True quantum discord, 
as a measure of quantum correlations must decrease with 
increasing temperature. 

It is interesting to notice that our treatment enables 
us to obtain many exact results exploring extensively 
the non-perturbative, non-Ohmic and non-Markovian 
regimes. In fact, analytic results for entanglement in the 
asymptotic regime can be obtained using only the long 
time values of the coefficients that appear in the exact 
master equation. Such values are known for a variety of 
physically relevant cases (see Ref. 



12j for a collection 



of such results). Numerical evaluation is needed only to 
take care of intermediate time regimes (until the depen- 
dence on the initial state is washed out by the evolution). 
This treatment is in contrast with previous works, such 
as the one reported in Ref. JLJ] that make extensive use 
of numerical solution. In particular, using our treatment 
it is rather simple to relax other previous assumptions 
(such as the restriction to low coupling or the limitation 
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to short time in the non-Markovian case) such as the ones 
usedinRef. [3]. 
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